home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / rsi_gammai.pro < prev    next >
Text File  |  1997-07-08  |  1KB  |  80 lines
  1. ; $Id: rsi_gammai.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6. ; modification: 09/21/92 - changed the number of iterations in
  7. ;               g_series (jiy)
  8.  
  9.  
  10. pro g_series,result,x,a
  11. ;evaluates incomplete gamma function with series representation.
  12. glog =lngamma(a)
  13. nsample=long(x/50.) > 1000l
  14. resarray = 1.0/(findgen(nsample) + a)
  15. resarray(1:*) = x*resarray(1:*)
  16. sum =1.0/a
  17. for i =1 ,nsample-1 DO BEGIN
  18. resarray(0) = resarray(0) * resarray(i)
  19. sum = sum + resarray(0)
  20. if (abs(resarray(0)) LT abs(sum)*3.0e-7) THEN BEGIN
  21.   result = sum * exp(-x + a * alog(x) - glog)
  22.   return
  23. ENDIF
  24. ENDFOR
  25.  
  26. result = -1
  27. return
  28. END
  29.  
  30.  
  31. pro g_fract, result,x,a
  32.  glog = lngamma(a)
  33.  gd =0.0 & fc =1.0 & b1= 1.0
  34.  bo = 0.0 & ao =1.0
  35.  
  36.  a1 = x
  37.  for n=1,100  DO BEGIN
  38.     an = float(n)
  39.     ana = an -a
  40.     ao = (a1 +ao * ana) * fc
  41.     bo = (b1 + bo *ana) * fc
  42.     anf = an * fc
  43.     a1 = x * ao + anf * a1
  44.     b1 = x * bo + anf * b1
  45.    if a1 THEN BEGIN
  46.      fc = 1.0/a1
  47.      g = b1 * fc
  48.      if abs((g-gd)/g) LT 3.0e-7 THEN BEGIN
  49.        result =  exp(-x + a* alog(x) - glog) * g
  50.        return
  51.      ENDIF
  52.     gd = g
  53.    ENDIF
  54.  ENDFOR
  55.  result =-1
  56. RETURN
  57. END
  58.  
  59.  
  60.  
  61. function rsi_gammai,a,x
  62.  
  63. if  x LT a+1.0 THEN BEGIN
  64.    g_series,result,x,a
  65.    return, result
  66. ENDIF ELSE BEGIN
  67.      g_fract,result,x,a
  68.      if result ne -1 then return, 1.0 - result  $
  69.     else return, -1
  70. ENDELSE
  71.  
  72. END
  73.  
  74.  
  75.   
  76.      
  77.  
  78.  
  79.  
  80.